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The Willmore flow is well known problem from the differential geometry. It min- 
imizes the Willmore functional defined as integral of the mean-curvature square 
over given manifold. For the graph formulation, we derive modification of the 
Willmore flow with anisotropic mean curvature. We define the weak solution 
and we prove an energy equality. We approximate the solution numerically by 
the complementary finite volume method. To show the stability, we re-formulate 
the resulting scheme in terms of the finite difference method. By using simple 
framework of FDM we show discrete version of the energy equality. The time 
discretization is done by the method of lines and the resulting system of ODEs 
is solved by the Runge-Kutta-Merson solver with adaptive integration step. We 
also show experimental order of convergence as well as results of the numerical 
experiments, both for several different anisotropics. 
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1. Introduction 

This article extends the isotropic Willmore flow defined by T. J. Willmore in 
[TT] . It is a minimizer of the Willmore functional defined as 



where T is given manifold smooth enough so that the mean curvature H can be 
evaluated almost everywhere on F. Prescribing the following normal velocity 
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we generate class of manifolds F (t) so that F (t) minimizes the Willmore func- 
tional as if: goes to infinity. One recognizes several formulation of the Willmore 
depending of the form in which F {t) is expressed - the graph formulation [U [HI [12] , 
the level- set formulation [51 [I] or the parametric formulation [liilj. Also approx- 
imation by the phase-field model exists [6] . In this article we will study only the 
graph formulation. Applications of the Willmore flow can be found in biology [3] 
or in image processing in image inpainting [2|. Especially in the later one, the 
anisotropic model can be useful. 



2. Problem formulation 

We assume having a manifold Fq described as a graph of function uq of two 
variables: 

Fo = {[x,Mo(x)]|xG^]cR'}, (1) 

where Q = (0, Li) x (0, L2) is an open rectangle and we will denote dQ its bound- 
ary. We assume having convex function 7 : ]R""'"^\{0} — )■ Rq , 7 = 7(^1, ■ ■ ■ Pn, — 1) 
which is positive 1-homogeneous i.e. 7 (Ap) = |A|7(p). We will call 7 surface 
energy density and denote 

where dp_'j stands for partial derivatives of 7 w.r.t. variable pi. We define 
anisotropic mean curvature of a manifold Fq smooth enough induced by energy 
density function 7 as 

H,{uo) = V-{Vp^{Vuo,-l)). (2) 
and the anisotropic Willmore functional as 

{uo) = 1 [ H'^ (uo) Q (uo) dx, (3) 



for Q (u) = yl + |Vm|^. We aim to find a function u* minimizing (3). The 
Euler-Lagrange equation for this functional takes the following form 



V.(^E,(«-)Vu,(«-)-2^V..j=0. 



(4) 



where we denoted 



(u) := Q (u) (u) , 
E^{u) := 9p^9p^.7(Vu,-l) = (Vp®Vp)7(Vu,-l). 
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In the rest of the text we will not emphasize explicitly the dependence of Q, H^, 
and on u. Multiplying ^ by the test function Lp G C°° and integrating 
over Q we get 



1 wl 

V ■ ( E^Vw^ - 2 0^^" ) 



E^Vw^ ■ + T^TjJVu ■ u^dS + / V ■ (E^Vw^) - - V ■ t^Vu ^dx. 
an Jn ^ \Q / 

The boundary integral over dQ vanishes if we set 

d„u = on dil, (5) 
E^Vw., -u = OondQ. (6) 

These relations define the Neumann boundary conditions. The Dirichlet bound- 
ary conditions 

u = (?! on dfl, (7) 
= g2 on dfl, (8) 

may be obtained in the same way but taking ip G (fl). In practice, the L2- 
gradient flow of ^ is solve rather than Q - see jH |5] . The following definition 
involves parabolic partial differential equation with unknown function u which is 
now also dependent on artificial time parameter t i.e. u = u{t,x.). In the same 
sense as in ([T|, a moving manifold T (t) is obtained. 

Definition 2.1. Let Q be a domain in R^. The anisotropic Willmore flow 
of graphs with the Dirichlet boundary conditions and the initial condition 
Uq is a fourth order parabolic problem given by 

f \ 

dtu = -QV ■ Ie^Vw^- --^Vuj on{0,T)xn, (9) 

= QH^ on (0, T) x fi, (10) 

U \t=0 = Uq on fi, (11) 

u = gi, = g2 on dfl. (12) 

The anisotropic Willmore flow of graphs with the Neumann boundary 
conditions and the initial condition uq is a fourth order parabolic problem given 
by Q)-([7ip and 

d^u = 0, E^Vw^ ■u = on dfl. (13) 
We also define the weak formulation: 
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Definition 2.2. Let Q be a domain in R^. The weak solution of anisotropic 
Willmore fiow of graphs with the Dirichlet boundary conditions 



u 



gi on dQ, 
g2 on dQ, 



is a couple u, : (0, T) — )■ (Q) which for each test function V9, ^ G Hq (Q) 
and a.e in (0,T) satisfies, 



Jn Q Jn 

i = - / Vp7-VCdx. 
Jn 



Q 



with the initial condition 



1 w"^ 

-^Vm ■ V^dx a.e. in (0, T) (14) 

(15) 



U 4=0= ^^0- 



(16) 



The weak solution of anisotropic Willmore fiow of graphs with homo- 
geneous Neumann boundary conditions 



on dVl, 
on (9f2, 



is a couple u, w : (0, T) — )• (Q) which for each test function (p,^ E {^) and 
a.e. in (0,T) satisfies (14)-(i5) and the initial condition (16). 

For the proof of the numerical stability we will need the following theorem: 



Theorem 2.3. For the solution u,w^ of (14)-(16) with the zero Dirichlet bound 

ergy e 

(dtu)" 



ary conditions the following energy equality holds: 

\2 



n 



1 d 

Q "^"+2dt 



n 



H^Qdx = 0. 



(17) 



Proof. We differentiate (15) with respect to t 



f ^^dx- / ^M|^dx+ / E,VdtU-V^ = for aU e G i/o' (^) 
Jn Q Jn Q Jn 



which follows from 
d 



dt 



Vp7(V«,-l)-Ve = 



i=l 
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Substituting = d^u in (14) and ^ = in (18) we have 

\2 



I (^^dx- f (E^Vw^) ■ V^iMclx + [ I'^Vu-Vdtudx = 0, (19) 
Jn Q Jn Jn^Q 



[ *^dx- f ^^dx+ / (E,va,.).v», = 

Jn W Jn ^ Jn 



(20) 



The sum of (19) and (20) gives 



^ Q + Q g2 +2Q3V^ Vd.^dx-O. 



Since Vc^^m ■ Vm = d^QQ, ([T]) turns to 

(9tM)^ Stw^w^ 1 w^^tg 



^n Q Q 2 Q- 

which is indeed what we wanted to show because 
1 d 



dx = 0, 



1 d _ dtw^w^ 1 w'^dtQ 
2 dt"^"^ 2 dt Q" ~ ~~Q 2 g2 • 



We will demonstrate the anisotropic Willmore flow of graphs on the following 
(an)isotropies: 



liso (P) 

7g (p, -1) 

labs (P) 



i + IpI, 

'1 + P^Gp, 



-f^^ + Cabs ^ -Pf 



(21) 
(22) 

(23) 



where G is symmetric positive definite matrix G G R^'^ and we use notation 
P = (p, —1). In fact, jiso represents the isotropic problem. The Wulf shapes W 
of given anisotropics defined by 



W= fl {xGR"|(x,q)<7(q)}, 
|q|=i 

are depicted on the Figure [TJ Simple calculations show that 



(24) 



H. 



H, 



7g 



1 + |Vm 
/GVm 

V 7g 



E 





1 




E^. 

/ISO 




(- 


~Q 


^( 


G- 


GVm 


7c; \ 


7C! 



~Q 



GVm 

7(G 



(25) 
(26) 



Figure 1: The Wulf shapes of 7^50, and jabs from left to right. 



It is difficult to express -^^7^^^ and E^^^^^ in some compact form and so we only 
show partial derivatives of •jabs with respect to pi and pj for i,j = 1, 2. 

3 

E^absPi , Pi p . 

/ ^ T 1 3 1 ^ = 1'2, 

i=l V -P,^ + Cafes Z]fc=l -Pfc V Pi + ^abs Zl.,=l -Pf 
3 

1abs,piPi — ^ ] ^ ^ 

1 »2 

+ ^ for i = 1, 2, 



Cabs ^absPi 



pl + ^absY.UP! {pl + eabsEUPl 

82 2 

_ ^absPiPj \ ^ ^absPiPj 

labs,pipj n 2-^ , q X 3 • 

^=1 {PI + eabs Ell PI) ' ^=1 {PI + eabs Ell ^,^) ^ 

3. Space discretization 

3.1. The complementary-finite volume method 

In [9] we applied the complementary-finite volume method for the space dis- 
cretization of the anisotropic surface-diffusion flow. In the same manner we dis- 



cretize even the problem from the Definition 2.1 Let /ii, /12 be space steps such 
that hi = and /i2 = ^ for some Ni, N2 G N"*". We define a numerical grid, 
its closure and its boundary as 

uJh = {{ihujh2)\t = l---Ni-l,j = l---N2-l}, (27) 
uJh = {(z/ii,j/i2) M = 0---iVi,j = 0---iV2}, (28) 
duh = uJj;\ojh- (29) 



For u E C (f2j we define its piecewise constant approximation on w as a grid 
function defined as {ihi,jh2) := m^- := u {ihi,jh2). We also define a dual 
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Figure 2: (a) Dual mesh for the complementary finite volumes method - dots denote ojh (27 1, 
circles denote dijJh ( 29 ) and solid lines stand for Vh ( 30 1 . (b) Notation for v^j j-j . 



mesh Vh as (see the Figure [2] a)) 



2) ^' V 2 

, = i...iv,-i,j = i...iV2-n. 



For < z < A^^i, < j < N2, i and j fixed, consider a finite volume Vij 
dual mesh Vh, denote its interior as Qij, its boundary as Tij and let /i (f^ij) 



(30) 



of the 



volume of ilij. We also denote all the neighboring volumes of the volume Vij as 
Afij. For all finite volumes Vij of the dual mesh Vh, 



the boundary Tij consists of 
four linear segments. We denote them as Fj^j-. It means that Fj^jj is a boundary 
of the finite volume Vij between nodes (i,j) and [i,])- By l^jj-j we denote the 
length of this part of Tij. We approximate the partial derivatives of u on the 
boundary of the finite volume Vij as (here resp. stands for the finite 
difference w.r.t. xi resp. X2): 











hi 






+1 ~ 




h2 



8^ 11^ 



Ui-lj 



hi 



ho 



and 



where we denote 

.h 



4 



4 

The quantity Q on the boundary of Fjj is approximated as 



Qij,i+lj ~ ^1 + + (^X2^ii,i+lj) ' 



Qij,ij+1 ~ \/l + + (^xa'^^ij'.ji+l) 



Qij,i-lj — \/l + + {9x2'^ij,i-lj) } 

and Q inside the finite volume 

cLS 

Integrating (|2]) over a finite volume iljj and applying the Stokes formula we obtain 



/" H^dx= I V-(Vp7)dx= I Vp7-z/dS. 



(31) 



where u denotes the outer unit normal vector to the finite volume boundary Tij. 
We approximate the term on the left as 



h 



(32) 



and the term on the right as 

I Vp7-i^dS= I Vp7-z/dS (33) 



For the inner finite volume Vij G Vh-, there are four different neighbors v-q G Mij 
(see the Figure 30). All the boundaries Fj^jj are linear segments and so v = Uij-fj 
is constant there. Moreover we assume that Vp7 is constant along Fj^jj too. It 
gives 

/ Vp7-z/dS?^i Y k^Y^^3J3■^i3^ (34) 
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where 



Vp7iiJi = {dp.lijjj, ^p,^^JJ]f = (5pi7 (VM-,-ij, -1) , dp.'y [Vu^^jj, -l)f. (35) 

(36) 



Putting ([32|) and (34) together we get 



For the regular dual mesh (30) we have 



= {h2^ ■ (1, 0)^ + /iiVp7ij- ij+i ■ (0, 1)^ 

+ /''2Vp7ij- ■ (-1, 0)^ + /iiVp7ij-ij_i ■ (0, -1)^ ^ 



( 9pi'Jij,i+lj 9pi1ij,i-l,j _j_ 9p2'yij,ij+^ ^P2lij,' 
hi h2 

and the approximation of = QH^ on the finite volume Qij is 

,h /^h Tjh 



W. 



We set 



w: 



h 



7:«J:«- 



2 
1 



Integrating (|9| over fijj and applying the Stokes theorem we get 



—d^dx 

Q 



1 wi^ 

KryVw-yU — - —^dj,udS 



(37) 

(38) 

(39) 
(40) 

(41) 



where u is the unit outer normal of the boundary Tij. The left hand side is 
approximated as follows: 



1 /i (%) d ^ Id;, 



(42) 

where we assume that u^j and Qf^- are constant on the element Vij. For the integral 
on the right hand side of (41) we have 

1 



with the usual notation. Putting (42) and (43) together gives 



Q 



v 



1 1 w'^ 



3 ^ 



In the terms of the regular dual mesh (30) we get 



hih2 



EV-,+i,V<,,,+i,.(1,0)^ 



V<,+i,.(l,0)' 



+ /i2 E: 



-1,0) 



T 



1 W. 



(Qi'j.j-ii) 



3 



-i,oy 



+ ^1 I E;,ij-ij_iVw;^^.,^, 



where for 



we have 



7>«i:«+ii 



7>11>U>U 7>12,ij,ij 

E ni • • -- E 



r 



(44) 



•^pi^piT (Vnf^-i+y, -1) , E^^-i^.+i = dp,dp,-f (Vu(;.ij.+i, -1) , 



The Neumann boundary condition di,u = on dQ takes the following discrete 
form 



if i = 1 then u 



(-1,0) = 


1 


i-h - 


<j) = 0, 


= (1,0) = 


1 


«. 


- ^JVi-l,i) 


(0,-1) = 


1 


i<l - 


uto) = 0, 


= (0,1) = 


1 







(45) 
(46) 
(47) 
(48) 
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and from (13) we get 



(49) 



Ni-lj,Nij 



j,Ni-lj,Nij 



^'y,12,Ni-lj,Nijdx2l^^,Ni-lj,Nij 



+ 



0, (50) 

E'7,21,jl,iO<9a.iW^ jl jO + E;^,22,il,i0f^a;2^7,jl,i0 ~ 0' 

(51) 



if i = 1 then z/ = (-1,0) 
a i = Ni-l then = (1, 0) 

if j = 1 then z/= (0,-1) 
if j = A^2 - 1 then u = (0, 1) 

3.2. The finite difference method 

We will prove the discrete version of the Theorem 2.3 in very similar manner 
as in [S]. For this purpose we need to re-formulate the complementary finite 
volume scheme by the finite difference method. We complete the numerical grid 
with virtual nodes representing the finite volumes edges and we denote the new 
grid by Wh- We index the nodes of Wh by indices i ^ \ and j ± ^. We define 
'^\±-j±- ~ \ (""iiijii + ""ij) ^"^^ denote by the uppercase a grid function with 
doubled indices i.e. 



0. 



(52) 



U'^i = u\i for A;,/ GO 

2 2 

We introduce the finite differences as follows: 



2N. 



(53) 



U 



f.,ki 



Tjh 



Tjh 



1 
2 



Tjh 


1,1 ^kl 


Tjh 

5 '-'b.,kl 

Tjh 

1 .b,kl 


jjh _ jjh 


Tjh 

^k,l 


h 

+1 '-^kl 


h 

jjh __ jjh 




h 


h 


■h 

f.,kl 


+ Kki) 


Tjh 

1 ^ .c,kl 





To approximate the gradient of u we define 

for i,jel---iV-l 



V 



l,'^c.,2j,2j> '-^.c,2j,2j; 



Note that since we define V^u^a in terms of Ub we can also write V^u^, i . , i for 



j = 1 ■ ■ ■ — 1. The terms in (|9|)-(|10|) are approximated as follows 



Tj-h 



C,r,G{-l,l) 

= ^7 («a 
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and the scheme has a form 



w. 



dt 

h 



(54) 
(55) 



Approximation of the Dirichlet boundary conditions is straightforward and the 
Neumann boundary conditions are handled as in (45)-(52). In what follows we 
define several discrete scalar products and we prove the discrete version of the 
Green formula. Assume having /, : cJ^ — )• R, / : aJ/i — )■ and the related 
functions F, G : tUh — )- R, -F : zuh — ?• R^ given by (53) we define 





4 ^ 

k=p,l=q 











2N-l,2N-l 



-' + [F,G,]Z^'^-\ 
'' + [F,G,]Z-'''^, 
{F,V'G)^ = {F\G.l+{F\G..)^ 

For the purpose of analysis, we will need the following grid version of the Green 
formula 



Lemma 3.1. Let u^,v'^ : -> M and U^, : Wh ^ ^ related to u,v by ( |5g| j. 

Then the Green formula is valid: 



+ 



+ 



- (n\vV) 



''2N-1 



2 



1=1 

IN~1 



k=l 



Proof. It is easy to see that for fixed k,l = 0, ■ ■ ■ 2A^ the following relations 
holds. 



\ri l^/^l^^"^'' _L h. (Tjh yh _ TThyh\ 



1 fc,2Af-l 



l^'^b.lii +9 {'-'2N,l^2N,l '-'iroi) 



k,2N-l h2 ^ h yh 

^ O \'~^ k,2N-l^ k,2N 



U^Vko) 



L ' -^Jfel ^ 9 \^k,2N^k,2N ^kl^ko) 



(56) 
(57) 
(58) 
(59) 
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Putting together (56) with (57) resp. (58) with (59) and summing over / 



1, ■ ■ ■ 2N — 1 resp. over = 1, ■ ■ ■ 2N — 1 we obtain 



(^in. = - 

which gives 



hi 



2N-1 



1=1 

2N-1 



k=l 



+ 



+ 



hi 

2 

h2 



(f/,VV)^ 

''2N-1 



1=1 

''2N-1 



E [K2N-1 + K2n) VI2N - Ho + ut) v;y 



k=l 



The following form of the previous Lemma 3.1 will be more convenient for us: 



Lemma 3.2. Let p , u , f : w/i — )■ M and v \duj= 0. Then 



(60) 



Theorem 3.3. For solution of (54)-(55) such that \du)h= and = \q^^ 
the following energy equality holds 



Proof. We start with the equation for w^^^^- (55), divide by Q^-, multiply by 
vanishing on duh and sum over uj. 



w: 
Q 



l^] =(V'^-(7.(W,-l)),0; 



The Green theorem (60) gives 

,h 



Q 



(v^e,Vp (vV,-i)). 



(61) 
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Taking the right hand side of (54), multiplying by a test function Lp vanishing at 
dujh and applying the Green theorem (60) we obtain 

,2 



h^3 



h^3 



(62) 



Differentiating (|61|) with respect to t we obtain 

,h 



d I w. 
Jt 



d 



U] +^(V,.7(W,-1),V'^0, 



I ( ^, e 1 + ((V, ® V,) 7 ( w, -1) vM, v^o, 



w. 



After substituting ,^ = iw^* we obtain 



w. 



Q 



l,t h 



(63) 



and a substitution (y9 = in (62) gives 



(Q 



/^^3 



0. 



(64) 



Substituting ( |63[ ) to ( |64[ ) (term E V w ) we have 



h\2 



1 



Q 



h ' "^7 



0. 



(65) 



We remind that V'^u^ ■ V^mJ* = ' which gives 



K) 



+ 



Q 



h ' "^7 



.,K) ) + 



0. 



(66) 
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It is equivalent to 



for 



K) 



Sh 



Finally from (67) we have 



{Q 




h \ ^ / h \ ^ 



h \ ^ / /i \ ^ 




t,iO 



+ 



(67) 



The term S*/! is vanishing on dco because of the zero Dirichlet boundary conditions. 



4. Time discretization 



The semi-discrete problem (44) is discretize in time by the method of lines 



in the same way as in [TO]. Here, the resulting system of ODEs is solved by the 
Runge-Kutta-Merson solver with adaptive integration time step implemented in 
CUBA to run on GPU. 



5. Experimental order of convergence 

To our best knowledge, there is not known any analytical solution for the 
Willmore flow of graphs even in the isotropic case. To be able to measure the ex- 
perimental order of convergence we will solve a modified problem with additional 
forcing term F having the following form: 



dtu 



-gv ■ ( E^vw^ 



2giVMj+F(M) on (0,T)xfi, (68) 

- QH^i on (0, T) X Vt, (69) 
u \t=Q = uq on Q, (70) 



u 



9, 



on dn. 



(71) 



Having a function (x, t), we set 

1 



Fw{C) 
(0 



Q(C)V 



P(C)V^, (0 



,Q(C)" 

Q{C)H^{0 onfix [0,r]. 



i ^^(C) , 

(0 



VC ]+dtCon^lx [0,T] 
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As an analytical solution C(x, t) of (68)-(71) we chose the following function 

ax,y,t) :=cos(7rt)-^(x"-r")(t/"-r")exp(-a(x2 + y2)) on Q x [0,T], 

(72) 

for Q = [— r, r]^. For given T, we evaluate the errors in the norms of the spaces 
Li {fl; [0,T]), L2 (fi; [0,T]) and L^o [0,T]) resp. their approximations 

Ni,N2,M 

ll^'-^'^(0||'';(^,;[o,T]) ■■= E r\u1^ikr)-C{-r + th,-r + jh,kT)\h', (73) 

i,j,k=0 



\U 



h''-'Ph(0\\LMo,T]) •= ^J^^^^Ju1^{kr)-C{~r + th,-r+jh,kr)\, 

fc=0,-- ,M 

for r = T/M. We would like to emphasize that r does not correspond with the 
time step of the solver. The experimental order of convergence is evaluated as 
follows - for two approximations u'^^ and u^^ obtained by the discretization with 
the space steps hi and h2 we compute the approximation errors Errh^ and Errh2 
in one norm of (73)-(75) and then we set 

Tpnr^fT? TP ^ log {Err h J Err h^) 

E0C{Errh„Errh2) := . . . • 76 

log(/ii//i2) 



The results are presented in the Tables [lj-[8j The solution of ([68j)-(j71j) was ap- 
proximated on the domain Q = [—4, 4f and the time interval [0, 0.1]. The Tables 
[l]-|2]show EOC for u and w with the isotropic setting. Both unknown functions 
are approximated with the second order. The Tables [s] - 14] show anisotropy (22) 
where 

o^(^;). ,rr) 

Even in this case we obtained EOC equals to 2. As we can see in the Tables [5]- 
|6} if we set 

' 2 1 



G = ^ ^ J , (78) 

mixed partial derivatives are involved and EOC drops drops to 1 for both u and 
w^. The last test is with the anisotropy given by (23) and we set eats = 0.1. 
Results are presented in the Tables [7] and [8j The non-linearity is much stronger 
here and it is not clear whether EOC is 1 or 2 from the Tables [7] and HI It would 
require computation on finer meshes which we were not able to do because of 
very long time needed for such simulation. 
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0.5 

0.25 

0.125 

0.0625 

0.03125 



Err. 

4.37- 10-3 
9.96 ■ 10-^ 
2.49 ■ 10"^ 
6.25 ■ 10-5 



EQC 

4.5 
2.13 
1.99 
1.99 



L2 (fi) 

Err. 

TL2 

4.39 ■ 10-3 
7.55 ■ 10-^ 
1.89 ■ 10-^ 
4.74 ■ 10-5 



EQC 

5.52 
2.54 
1.99 
1.00 



^00 m 

Err. 

X05 

2.82 • 10-2 
1.68 • 10-3 
4.22 • 10"^ 
1.05-10-^ 



Table 1: EOC of the approximation of the function u with the anisotropy defined by (21) 



0.5 

0.25 

0.125 

0.0625 

0.03125 



Err. 

TU8 

3.19 ■ 10-2 
7.02 ■ 10-3 
1.76 ■ 10-3 
4.40 ■ 10"^ 



EOC 

3.92 
2.18 
1.99 
2 



L2 iS^) 

Err. 

7L83 

4.89 ■ 10-2 
7.79 ■ 10-3 
1.91 ■ 10-3 
4.75 ■ 10"^ 



EQC 

4.09 
2.65 
2.02 
2 



Loo 

Err. 
T65 
0.53 

8.81 • 10-2 

2.1 • 10-2 

5.2 • 10-3 



Table 2: EOC of the approximation of the function with the anisotropy defined by (21 ) 



Meshes 



h 



Err. 



EOC 



L2 (O) 



Err. 



EOC 



Loo m 



Err. EOC 



16 
32 
64 

128 
256 



0.5 

0.25 

0.125 



0.7 

0.53 

0.085 



0.0625 0.0059 
0.03125 0.0015 



0.4 
2.6 
3.8 
2 



1.2 

0.93 

0.24 

0.013 

0.0035 



0.37 
2 

4.2 
1.9 



5.5 

4.9 

1.5 

0.089 

0.024 



0.17 
1.8 

4 

1.9 



Table 3: EOC of the approximation of the function u with the anisotropy defined by ( 22 ) and 

(TttI). 
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Meshes 


h 




(^]) 


L2 m 


Loo 




Err. 


EOC 


Err. 


EOC 


Err. 


EOC 


16 


0.5 


3.5 




7.3 




41 




32 


0.25 


3.2 


0.15 


7.3 





44 


-0.1 


64 


0.125 


0.58 


2.4 


1.6 


2.1 


19 


1.2 


128 


0.0625 


0.04 


3.9 


0.085 


4.3 


0.5 


5.2 


256 


0.03125 


0.01 


1.9 


0.023 


1.9 


0.14 


1.9 



Table 4: EOC of the approximation of the function with the anisotropy defined by (22) and 



Meshes 


h 


LiiSl) 


L2 


m 


Loo 




Err. 


EOC 


Err. 


EOC 


Err. 


EOC 


16 


0.5 


2.1 




4 
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32 


0.25 


0.46 


2.2 


0.64 


2.6 


2.9 


2.7 


64 


0.125 


0.15 


1.6 


0.34 


0.92 


1.9 


0.59 


128 


0.0625 


0.06 


1.3 


0.16 


1.1 


1.1 


0.74 


256 


0.03125 


0.025 


1.2 


0.08 


0.99 


0.72 


0.66 



Table 5: EOC of the approximation of the function u with the anisotropy defined by ( 22 ) and 
pi. 



Meshes 


h 


Li 




L2 


(^]) 


Loo 








Err. 


EOC 


Err. 


EOC 


Err. 


EOC 


16 


0.5 


12 




27 




120 




32 


0.25 


3.3 


1.9 


5 


2.4 


30 


2 


64 


0.125 


1.2 


1.5 


2.7 


0.9 


25 


0.26 


128 


0.0625 


0.53 


1.1 


1.3 


0.99 


13 


0.92 


256 


0.03125 


0.25 


1.1 


0.76 


0.82 


8.3 


0.68 



Table 6: EOC of the approximation of the function with the anisotropy defined by (22) and 
(l78|. 



18 



Meshes 


h 


Li in) 


L2 




Loo 




Err. 


EOC 


Err. 


EOC 


Err. 


EOC 


16 


0.5 


2.9 




3.2 




9 




32 


0.25 


1.8 


0.72 


1.7 


0.9 


4.4 


1 


64 


0.125 


0.46 


1.9 


1.5 


0.17 


15 


-1.8 


128 


0.0625 


0.078 


2.6 


0.11 


3.7 


0.55 


4.8 


256 


0.03125 


0.027 


1.47 


0.045 


1.28 


0.18 


1.6 



Table 7: EOC of the approximation of the function u with the anisotropy defined by (23) and 

Cafes = 1- 



Meshes 


h 


L 


i(^) 


L2 




Loo 




Err. 


EOC 


Err. 


EOC 


Err. 


EOC 


16 


0.5 


14 




14 




41 




32 


0.25 


8.6 


0.69 


8.3 


0.74 


29 


0.49 


64 


0.125 


9 


-0.067 


36 


-2.1 


460 


-4 


128 


0.0625 


0.69 


3.7 


1.1 


5 


5.3 


6.4 


256 


0.03125 


0.29 


1.25 


0.49 


1.16 


1.49 


1.8 



Table 8: EOC of the approximation of the function with the anisotropy defined by (23) and 

Cafes = 1- 
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6. Numerical experiments 



On the Figures |3]-[5] we show results of qualitative analysis. The initial con- 
dition is uq = sin ^tta/x^ + j on 1] = [—2,2]^ and we set the Neumann 
boundary conditions d^ip = E^Vw-y ■ u = on dQ. The computational domain 



is covered by 100 x 100 meshes. The Figure [3] depicts result obtained with the 
anisotropy given by (22) and 



8 
1 



(79) 



on the time interval [0,0.001] . The Figure |4] show the same anisotropy (22) but 
with 

« - ( 8 10 ) («°) 

on the time interval [0, 1.024 ■ 10"'^]. And finally the Figure [s] reveals the behavior 
of the anisotropy given by (23) with eabs = 0.001 on the time interval [0, 0.006]. 



7. Conclusion 



We have presented new mathematical formulation for the anisotropic Will- 
more fiow of graphs. We have proved energy equality for the new problem. To 
approximate the solution numerically we showed complementary finite volume 
based scheme. We also show how to re-formulate the numerical scheme by the 
finite difference method. This approach allowed us to prove stability of the nu- 
merical scheme. The computational part of the article contains experimental 
order of convergence evaluated on both studied anisotropies. We also showed few 
qualitative results. 
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